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Air  Force  by  providing  more  accurate  assessments  for  systems 
safety  purposes.  This  could  reduce  unnecessary  costs  in  ground 
operations  and  enable  an  increase  in  sortie  generation  rate.  The 
developed  system  would  be  applicable  to  other  military  services 
as  well  as  other  fuels,  such  as  missile  fuels,  and  could  be 
useful  for  system  safety  studies  for  commercial  operators. 
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1.0  INTRODUCTION 


The  loading  and  off-loading  of  fuel  from  aircraft  is  a  process  in 
which  a  significant  system  failure  can  occur.  Such  a  failure 
could  result  in  a  release  of  flammable  liquid  and  vapor.  The 
migration  of  the  potentially  hazardous  fuel  vapors  pose  a  threat 
to  the  safety  of  nearby  personnel  as  well  as  to  ground  support 
equipment  and  other  aircraft.  When  considering  safety  and 
economics  of  ground  support  operations,  the  accurate  prediction 
of  the  physical  behavior  of  the  fuel  releases  is  extremely 
important.  While  overly  conservative  and  simplistic  simulation 
models  could  be  used  to  provide  assessment  of  personnel  safety 
and  equipment  vulnerability  relative  to  the  fires  resulting  from 
fuel  vapors,  these  types  of  models  can  suggest  mitigations  that 
are  unreasonable  and  which  could  seriously  restrict  ground 
support  operations. 

Tracer  Technologies  understands  that  the  development  of  models 
which  can  accurately  predict  all  possible  accident  events  would 
be  extremely  expensive  and  unrealistic.  Private  industry  has 
invested  much  time  and  money  during  the  past  decade  toward  the 
development  of  dispersion  models/methodologies  which  can  predict 
vapor  migration/dispersion  from  various  spills  of  flammable 
materials.  Personnel  at  Tracer  Technologies  were  responsible  for 
the  development  of  many  of  these  models  and  calculational  methods 
which  have  been  successfully  employed  by  the  oil  and  gas 
industry.  The  purpose  of  this  study  was  to  modify  and  add  to 
these  models  in  order  to  derive  a  family  of  models  which  could 
be  used  routinely  by  the  U.S.  Air  Force  to  study  the  potential 
consequences  relative  to  leakages  and  spills  of  fuels.  The 
models  would  then  provide  input  to  system  safety  engineering 
analysis  and  promote  the  development  of  cost  effective  and 
optimal  material  handling  procedures  for  ground  support 
operations. 
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The  technical  objectives  of  the  program  as  outlined  in  the 
original  proposal  to  the  Air  Force  were: 

1  -  Define  the  types  and  sizes  of  fuel  spills  expected  at 
military  ground  support  operations. 

2  -  Define  the  range  of  conditions  under  which  the  fuel  spills 
occur  at  military  bases  and  ground  support  operations. 

3  -  Develop  or  modify  existing  Tracer  Technologies  numerical 
models  to  define  time  dependent  vapor  source  rates.  These  rates 
will  be  dependent  upon  fuel  characteristics,  storage  and  transfer 
methodologies,  and  environmental  conditions. 

4  -  Execute  the  models  in  an  analysis  designed  to  aid  in 
identifying  the  key  physical  parameters  which  govern  the  behavior 
of  various  types  of  spills  envisioned. 

5  -  Modify  and  develop  two  dimensional  dispersion  models  capable 
of  providing  realistic  results  for  small  to  medium  spills  that 
can  occur  at  Air  Force  ground  support  operations. 

SUMMARY  OF  PROGRAM  RESULTS 

The  five  objectives  noted  have  been  completed  with  particular 
emphasis  of  JP4  fuel.  This  material  is  the  most  volatile  of  the 
fuels  presently  in  use  in  military  aircraft  and  presents  the 
greatest  hazard.  This  emphasis  was  agreed  upon  during  initial 
technical  meetings  with  Air  Force  experts.  A  review  of  the 
possible  kinds  of  releases  that  could  occur  during  typical 
refueling  operations  indicated  that  all  of  the  releases  could  be 
fit  into  the  general  catagories  of  (1)  vent  spills  or  releases 
which  occur  during  refueling  as  gaseous  releases,  (2)  liquid 
pool  spills  which  would  occur  if  a  fuel  line  failed  and  liquid 
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material  was  spilled  on  the  ground,  and  (3)  pressurized  hose 
breaks  which  could  result  in  the  spraying  of  liquid  droplets  into 
the  atmosphere.  These  types  of  scenerios  were  then  considered 
as  the  generic  kinds  of  spills  which  would  be  used  in  the  model 
development  portion  of  the  contract. 

An  examination  of  the  literature  and  discussions  with  personnal 
at  both  Norton  AFB  and  Wright  Patterson  AFB  have  served  to  define 
probable  spill  types  and  sizes  for  this  study.  Ciccone  and  Graves 
(  1  ) ,  in  a  study  to  determine  the  methods  for  neutralization  of 
DOD  aircraft  fuel  spills,  determined  that  spills  could  be 
catagorized  as  small,  medium  or  large  as  follows;  small:  0-4 
gallons,  medium:  5-42  gallons,  and  large:  greater  that  42 
gallons.  Discussions  with  present  Air  Force  engineering 
personnel  indicate  that  present  operating  procedures  would  allow 
failures  of  approximately  this  size  given  a  specific  scenerio. 
For  example,  a  3"  diameter  fueling  hose  50  feet  long  (typical  of 
that  used)  contains  approximately  18.4  gallons.  Allowing  an 
operator  a  delay  time  of  1  to  5  seconds  to  respond  to  an 
emergency  and  cause  system  shut  down,  a  total  spill  of  26  gallons 
to  56  gallons  would  occur.  For  purposes  of  the  calculations 
performed  during  phase  1  of  this  study,  two  spill  sizes  were 
selected,  4  gallons  and  40  gallons.  The  smaller  spill  was  chosen 
since  it  represented  the  upper  bound  of  the  most  probable  spill 
grouping,  while  the  larger  spill  represented  the  a  median 
estimate  between  the  "medium”  and  "large"  spill  categories. 

The  spill  of  fuel  could  be  expected  to  occur  outside  or  within 
ground  shelters.  In  the  former,  all  meteorological  conditions 
must  be  taken  into  account.  In  the  latter,  low  wind  speed 
conditions  will  predominate.  Therefore,  during  this  study,  two 
dispersion  models,  each  representing  one  of  these  cases,  was 
developed  to  accurately  assess  Air  Force  needs.  Models  were 
developed  for  the  source  characterization  as  well  as  for  the 
transport  and  dilution  of  the  resulting  vapor  clouds  in  the 
atmosphere.  These  models  can  be  applied  to  specific  situations 


to  evaluate  the  potential  hazard  area  resulting  from  fuel 
spills . 

The  two  source  models  developed  were  (1)  a  liquid  pool  spill 
model  which  includes  time  dependent  heat  and  mass  transfer 
effects,  and  (2)  liquid  droplet  spray  model  which  includes  time 
dependent  particle  dynamics  as  well  as  heat  and  mass  transfer 
effects.  The  transport  and  dispersion  models  include  (1)  a 
steady  state  dispersion  model  describing  the  diffusion  of  vapor 
from  a  steady  source  under  moderate  wind  conditions  and  any 
atmospheric  stability,  and  (2)  a  time  dependent  vapor  dispersion 
model  specifically  for  low  speed  wind  conditions.  The  model 
theory  is  described  in  section  2.0  of  this  report  while  the 
actual  codes  are  discussed  in  Appendices  A  through  D.  Appendix  E 
contains  a  review  of  the  parameters  utilized  in  the  model 
calculations  presented  and  Appendix  F  provides  sample  outputs 
from  the  source  models  as  well  as  the  dispersion  models. 

The  models  were  developed  in  such  a  manner  that  data  input 
requirements  can  vary  from  very  simple  to  complex.  In  the  case 
of  minimum  initial  data,  only  the  spill  scenario,  the  spill  size, 
a  wind  speed  at  a  specific  height,  and  the  atmospheric  stability 
need  be  specified.  Should  it  be  available,  complete  atmospheric 
wind  speed  profiles,  spray  droplet  size  distributions,  and  source 
specifics  such  as  fuel  line  pressure,  hole  size  ,  etc.  may  be 
included  in  the  calculations. 

Tracer  Technologies  has  proposed  a  phase  II  of  the  program  which 
would  consist  of  both  continued  development  of  the  models, 
including  the  completion  of  a  modified  set  of  models  which  can  be 
programmed  on  hand  held  calculators  for  field  use,  and  an 
experimental  program  which  will  add  significant  credibility  to 
the  close  in  dispersion  coefficients  required  to  accurately 
assess  fuel  spills.  These  proposed  future  study  requirements 


2.0  TECHNICAL  DISCUSSION 
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Modeling  the  failure  of  an  aircraft  fuel  spill  is  done  in 
two  discrete  steps.  The  first  requires  the  determination  of  the 
source  rate  of  flammable  gases  into  the  atmosphere.  The  second 
step  is  the  prediction  of  the  dispersion  of  the  resulting 
flammable  vapors.  These  two  steps  will  yield  the  concentration 
of  the  fuel  as  functions  of  both  time  and  position.  One  can  then 
analyze  the  results  to  determine  the  extent  of  the  safety  hazard. 

In  order  to  calculate  the  source  rate,  Tracer  Technologies  has 
developed  two  models  which  account  for  tank  or  pressurized  line 
spills.  Pure  vapor  sources  are  directly  input  into  the  disper¬ 
sion  models.  Since  the  primary  hazard  associated  with  JP4  fuel 
is  a  fire,  the  models  need  to  deal  only  with  minimum  concentra¬ 
tions  on  the  order  of  1%.  The  lower  flammability  limit  for  JP4 
is  1.3%  (2).  Should  toxicity  become  a  major  concern  with  other 
types  of  fuels,  much  lower  concentrations  become  important  and 
changes  to  the  developed  models  would  be  required. 

The  source  models  require  a  description  of  the  physical 
parameters  and  conditions  associated  with  the  JP4  tank  or  fueling 
system.  For  tanks,  geometric  dimensions  such  as  tank  height, 
diameter,  and  any  confinement  devices  such  as  dikes  that  may  be 
in  the  immediate  area  are  of  importance.  For  fueling  systems, 
the  physical  parameters  are  required  as  well  as  line  pressures, 
temperatures,  and  delivery  rates  of  the  incoming  fuel.  Vent 
vapor  flows  require  only  the  flow  rate  and  exit  area  be 
described.  For  non-pressur ized  flows,  the  fuel  will  exit  as  a 
liquid  and  form  a  pool.  For  pressurized  flows,  the  fluid  may 
exit  as  a  gas,  a  liquid,  and  thus  predict  time  dependent  vapor 
flow.  Environmental  conditions  such  as  air  fuel,  and  ground 
temperatures  and  wind  speeds  are  also  required.  Should  they  not 
be  specified,  default  values  are  built  into  the  models.  The 
vapor  flow  rate,  given  as  a  function  of  time,  is  input  into  one 


of  the  two  dispersion  models  developed.  These  models,  developed 
for  low  or  moderate  wind  speeds,  require  a  description  of  the 
environmental  conditions.  Specifically,  a  wind  speed  and 
associated  height  and  a  description  of  the  atmospheric  stability 
is  required.  The  models  are  capable  of  taking  into  account 
complete  wind  profiles  should  that  information  be  available.  The 
downwind  dispersion  of  the  incoming  vapor  source  is  calculated  as 
a  function  of  both  time  and  position. 


2.1  Source  Model  Requirements 


The  vapor  source  models  require  that  the  evaporation  rate  of  JP4 
fuel  droplets  or  pools  be  described.  The  estimate  involves  a 
description  of  the  total  mass  transfer  process  from  the  bulk  of 
the  liquid  to  the  atmosphere,  the  heat  transfer  between  the 
surrounding  air,  the  ground,  and  the  droplet  or  pool,  and  the 
loss/gain  of  energy  from  the  liquid  to  the  atmosphere  through 
radiation  processes.  In  addition  to  the  evaporation  terms,  the 
dynamics  of  a  droplet  from  a  high  pressure  line  leak  must  be 
described,  as  well  as  the  size  and  thickness  of  any  liquid  pool 
forming.  The  description  of  the  processes  in  this  section  are 
organized  as  follows:  initially  the  heat  and  mass  transfer 
processes  are  discussed  since  they  are  common  to  both  models; 
secondly,  the  liquid  pool  model  is  be  discussed;  and  finally, 
the  vapor  source  associated  with  a  break  in  a  pressurized  fuel 
spill  is  described.  The  physical  parameters  that  describe  JP4 
are  important  to  all  of  these  calculations.  JP4  fuel  properties 
were  obtained  from  reference  (2),  "Handbook  of  Aviation  Fuel 
Properties".  The  data  was  complied  by  the  CRC  Aviation  Handbook 
Advisory  Group  from  the  latest  known  sources  on  each  particular 
subject.  Where  conflicts  arose  owing  to  discrepancies  in  source 
material,  they  were  resolved  by  decision  of  the  Group.  The  data 
was  published  in  1983. 


2.1.1  Evaporation  Rates  of  JP4  Fuel 

The  problem  considered  is  that  of  determining  the  rate  of 
evaporation  of  JP4  as  a  function  of  temperature,  wind  speed, 
atmospheric  conditions,  solar  radiation,  the  dimensions  of  the 
pool  or  droplet,  and  the  volatility  and  diffusion  characteristics 
of  the  fuel.  The  methodology  is  similar  to  that  utilized  by 
MacKay  and  Matsugu  (3).  This  estimate  requires  a  description  of 
the  total  mass  transfer  process  from  the  liquid  to  the 
atmosphere.  The  mass  transfer  process  can  include  the  liquid 
phase  mass  transfer  resistance,  which  controls  the  rate  of 
transfer  of  material  from  the  liquid  to  the  interface  and  depends 
on  the  dif f usivities  and  the  flow  conditions  in  the  liquid,  and 
the  transfer  from  the  interface  to  the  atmosphere.  In  this 
problem,  only  the  transfer  from  the  interface  to  the  atmosphere 
is  of  importance  since  liquid  concentration  gradients  are  assumed 
to  be  insignificant.  The  vapor  composition  is  assumed  to  equal 
that  of  the  liquid  and  the  liquid  composition  is  assumed  to 
remain  constant  with  time.  A  vapor  phase  mass  transfer 
coefficient  is  defined  in  equation  (1)  which  gives  the  mass 
flux  as  a  function  of  the  vapor  pressure  driving  force  which  is 
the  difference  between  the  JP4  vapor  pressure  at  the  surface  P 
and  the  vapor  pressure  in  the  atmosphere  P,,, 


(1)  N  =  Km(P-P.)/RTliq 


The  value  of  P  is  not  well  known.  Comparisons  of  the  theory 
with  droplet  evaporation  experiments  (ref  4)  were  made  and  a 
value  of  P  inferred  which  resulted  in  good  agreement  between 
the  two.  This  procedure  is  described  in  more  detail  in  Appendix  E 
of  this  report . 

K  is  a  function  of«the  transport  conditions  in  the  atmosphere 
m 

immediately  above  the  spill  or  surrounding  the  droplet. 
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2. 1.1.1  Heat  Transfer 

An  evaporating  drop  or  liquid  pool  generally  has  a  temperature 
close  to  that  of  the  atmosphere  and/or  ground.  At  high  rates  of 
evaporation,  JP4  may  be  considerably  cooler  due  to  evaporative 
enthalpy  losses.  It  may  warm,  however,  due  to  absorption  of 
solar  radiation,  or  heat  conduction  from  the  ground.  To  permit 
accurate  predictions  of  evaporation  rates,  these  effects  must  be 
quantified  and  the  liquid  temperature  related  to  local 
conditions.  The  liquid  temperature  is  determined  by  the  direct 
heat  transfer  from  the  air,  from  the  ground  if  considering  a 
pool,  the  incident  solar  radiation,  the  emitted  pool  radiation, 
the  rate  of  evaporation,  the  enthalpy  of  evaporation  and  the 
depth  and  temperature  history  of  the  pool  or  droplet.  For  the 
model  developed  here,  the  pool  or  droplet  is  considered  to  be  of 
uniform  temperature. 

2. 1.1. 2  Velocity  Characteristics  of  the  Lower  Atmosphere 

The  velocity  profile  above  the  JP4  pool  and  surrounding  the 
droplets  are  critical  in  controlling  the  eddy  diffusivity  and  the 
droplet  drag.  It  was  convenient  to  follow  Sutton  (5)  and  assume 
that  the  wind  velocity  profile  follows  a  power  law  as  indicated 
below  in  which  U  is  the  windspeed,  U1  is  the  windspeed  at  a 
specific  height  (1  meter)  and  Z  is  the  height. 


n 

<2>  U  =  U^2-11 


The  exponent  n  is  a  function  of  ground  roughness  and  temperature 
profile  in  the  atmosphere  (stability).  Typical  values  are  from 
0.25  to  1.00.  For  average  atmospheric  conditions  a  value  of  0.25 
for  n  is  reasonable  (5).  Sutton  has  shown  that  using  the  Von 
Karman  similarity  principle  and  this  assumed  velocity  profile, 
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the  turbulent  diffusivity  can  be  expressed  as  a  function  of  the 
velocity  and  the  mass  transfer  coefficient  can  be  deduced  as 
equation  1  where  C  is  a  constant  and  X  is  the  pool  or  droplet 
diameter.  C  is  taken  from  experimental  data. 


(3)  K  =CU°'78X"°*  11 

til 


_  2. 1.1. 3  Theory 


From  the  heat  and  mass  transfer  analogy  (6),  the  heat  and  mass 
transfer  coefficients  can  be  related  as  follows: 


(4)  (Stm)  (Sc)0’67  =  <sth)  (Pr)0’67 
thus 

Kh  '  VV^r'0'67 

The  rate  of  heat  transfer  from  the  atmosphere  is  K  (T  -  T,  )A 

h  A  liq 

cal/hr  where  T  and  T  are  the  atmospheric  and  liquid 

A  liq 

temperatures  and  A  is  the  surface  area.  The  heat  transfer  from 

the  ground  is  similiarly  expressed  as  K  (T  -T  )A  where  K 

c  G  liq  c 

is  the  heat  transfer  coefficient  from  the  ground  to  the  liquid, 
and  T^  is  the  ground  temperature. 

For  liquid  pools,  solar  radiation  is  calculated  as  the  incident 
direct  solar  radiation  and  scattered  sky  radiation  (L),  which 
reaches  the  earth’s  surface.  A  proportion  of  this  is  reflected, 
this  being  quantified  by  a  surface  albedo  "a".  The  heat  input 
from  solar  radiation  will  thus  be  (l-a)LA  cal/hr. 


loss  by  long  wave  radiation  from  the  pool  is  T  A 


where  e is  the  emissivity  of  the  pool  for  long  wavelength 
radiation  and  a  is  the  Stef an-Boltzmann  constant.  This  term  in 
the  equation  was  only  taken  into  account  for  liquid  pool  spills. 


Evaporative  cooling  will  cause  an  heat  loss  of 


(5)  K  AWAH  AP/RT 

n 


where  W  is  the  molecular  weight  of  the  JP4  fuel,  AH  the 
enthalpy  of  vaporization  and  AP  the  vapor  pressure  driving 
force . 

Combining  these  heat  flux  terms  with  the  mass  transfer  effects  of 
evaporating  JP4,  the  differential  equation  describing  the 
temperature  variation  of  the  perfectly  mixed  pool  with  time  can 
be  written. 


C6I  dTu, 


(l-a)L-eoT*  ,i+e,0TaWx,(pvC  „(SC/Pt>  — '  IT,-  1 


0.67 


p  1  tq  a 


-WAHV AP/RT  )  +  Kc<Tg-Tiiq) 


The  equation  is  integrated  numerically  to  determine  T^^as  a 
function  of  time.  Once  T  is  known,  the  rate  of  evaporation 
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and  cumulative  mass  losses  are  calculated  using  equation  3  and 
the  mass  conservation  equation. 

2.1.2  Liquid  Pool  Definition 

The  physical  parameters  which  describe  a  liquid  pool  were 
determined  using  data  from  Ciccone  and  Graves  (1).  In  their 
study  of  fuel  tank  leak  classification,  they  measured  the  surface 
area  of  spills  as  a  function  of  the  overall  spill  size.  An 
analysis  of  the  data  presented  lead  to  the  conclusion  that  the 
thickness  of  the  liquid  pool,  for  spill  sizes  considered  in  this 
study,  would  be  approximately  0.0005  meters.  This  was  treated  as 
a  constant  in  the  pool  source  rate  model  with  the  pool  area  then 
calculated  from  the  spill  volume.  Equation  5  was  integrated  and 
the  temperature  of  the  liquid  was  assumed  constant  through  out. 


2.1.3  Droplet  Dynamics 

JP4  fuel  emitted  from  a  broken  or  burst  fueling  line  exits  as 
liquid  droplets.  The  dynamics  of  the  droplets  are  necessary  to 
describe  the  vapor  source.  Figure  2.1  shows  the  forces  acting 
on  the  drop  as  it  passes  through  the  atmosphere.  The  equations 
describing  the  motion  are 
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(7)  Dy-W 
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dv 
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where  the  drag  force  (assumed  for  a  sphere)  is  expressed  as 

(8)  D  -  Cn  Ip  aV2A 
D  2  8 


( 


11 


Uh 


V 


vertical  velocity  v 

t 


horizontal  velocity  u 


weight 


Figure  2-1 

Dynamics  on  JP4  Fuel  Droplet 


and  can  be  positive  or  negative,  V  is  the  total  velocity  of  the 
droplet,  and  Ag  is  the  droplet  frontal  area.  The  weight  is 
expressed  as 


(8)  W  =  p , VOL  .  g  =  4*r3  ,  p.g 
1  drop’  3  drop  r 


The  equations  are  integrated  numerically  for  U  and  V,  with  the 
droplet  evaporation  accounted  for  using  equations  (1)  and  (6)  to 
calculate  a  new  droplet  radius  and  mass  evaporation. 

2.2  Modeling  Atmospheric  Dispersion 

With  the  source  rate  calculation  complete,  attention  is  turned  to 
dipsersion  of  the  vapors.  The  ambient  wind  will  begin  to  advect 
the  vapor  cloud  downwind.  Diffusional  forces  will  act  to 
disperse  the  vapors  outward  diluting  the  cloud  and  buoyancy 
forces,  if  present,  will  cause  the  cloud  to  rise  or  slump  adding 
to  the  dispersion.  All  of  these  factors  must  be  accounted  for  in 
order  to  accurately  predict  concentrations  at  locations  away  from 
the  failure  as  functions  of  time. 

Atmospheric  diffusional  forces  are,  as  one  would  expect,  strongly 
influenced  by  the  local  wind  field  and  the  stability  of  the 
atmosphere.  The  wind  guides  the  direction  of  the  cloud  movement, 
and  its  strength  influences  the  degree  of  air  entrainment  in  the 
vapor  cloud.  The  stability  of  the  atmosphere  is  the  tendency  of 
buoyancy  forces  to  enhance  or  oppose  vertical  motion.  The 
atmospheric  stability  is  a  function  of  many  factors,  but  one  of 
the  most  important  is  surface  or  ground  heating.  Hence,  the  type 


of  surface;  sand,  asphalt,  concrete,  etc.  and  degree  of  cloud 
cover  can  be  important  parameters. 

Wind  Conditions 


Wind  speeds  are  the  primary  factor  in  determing  the  transport  of 
the  vapor  cloud.  One  of  the  two  major  effects  of  the  wind  is  to 
advect  the  cloud  downwind  moving  the  hazardous  zone.  Higher 
winds  will  move  the  cloud  greater  distances  which  can  be  a  very 
negative  factor  from  a  safety  standpoint.  The  higher  winds  can, 
however,  have  a  beneficial  effect  in  that  higher  winds  generate 
higher  turbulent  levels  in  the  atmosphere.  The  greater  the 
turbulence,  the  greater  the  dispersion  diluting  the  cloud  down  to 
safe  concentration  levels. 

Specifying  wind  conditions  can  be  a  very  difficult  task.  Wind 
speed  and  direction  often  change  considerably  from  one  hour  to 
the  next.  The  effect  of  these  fluctuations  is  to  spray  the 
flammable  material  over  a  much  wider  area  than  would  be  predicted 
by  specifying  a  single  direction.  This  tends  to  make 
concentration  predictions  far  downwind  conservative  due  to  the 
increased  dilution. 


Atmospheric  Stability 


More  important  than  turbulent  dispersion  produced  by  eddies 
caused  by  the  wind  are  those  due  to  convection  from  density 
gradients.  This  convection  dispersion  is  strongly  influenced  by 
surface  heating  and  the  vertical  temperature  distribution  in  the 
atmosphere.  When  air  is  heated  from  below,  upward  moving 
convective  currents  develop.  These  are  caused  since  the  lower 
air  is  hotter  and  hence  more  buoyant  than  the  air  above.  When 
these  buoyancy  forces  enhance  the  verticle  mixing,  the  atmosphere 
is  described  as  unstable.  If  the  opposite  occurs,  the  colder 
air  tends  to  oppose  vertical  motion.  In  this  situation  the  air 
is  considered  stable. 


The  stability  of  the  atmosphere  can  be  characterized  by  examining 
the  rate  at  which  the  air  temperature  changes  as  one  moves  up  in 
altitude.  This  rate  of  temperature  change  is  known  as  the  lapse 
rate.  In  dry  air,  if  an  air  parcel  is  displaced  vertically  and 
is  always  in  equilibrium  with  its  surroundings  (hence  adiabatic), 
the  lapse  rate  is  approximately  1  C/100  m.  Lapse  rates  less 
than  this  tend  towards  instability  while  those  greater  tend  to  be 
more  stable.  Figure  2.2  shows  examples  of  these  graphically. 
The  stability  of  the  atmosphere  is  generally  broken  into  seven 
classes  based  on  the  lapse  rate.  These  are  specified  by  the 
letters  A  through  G  corresponding  to  extremely  unstable  to 
extremely  stable.  Table  2.1  shows  this  breakdown. 

Modeling 

The  traditional  method  for  modeling  the  dispersion  of  a 
substance  in  the  atmosphere  is  to  assume  diffusion  coefficients 
are  independent  of  position  and  that  the  diffusion  follows  Fick's 
Law.  By  further  assuming  that  the  wind  moves  at  an  average 
speed,  also  independent  of  position,  one  arrives  at  the  well- 
known  gaussian  equation.  This  predicts  a  binormal  concentration 
distribution  transverse  to  the  wind  as  one  moves  downwind  of  the 
source.  One  predicts  the  diffusion  coefficients  based  on  the 
stability  of  the  atmosphere. 

This  type  of  modeling  is  fairly  good  provided  the  source  in 
neutrally  buoyant  and  the  other  assumptions  are  met.  However, 
when  dealing  with  the  type  of  releases  that  occur  from  pipeline 
breaks  or  tank  failures  this  type  of  modeling  in  inappropriate. 
As  discussed  earlier,  the  wind  field  is  not  independent  of 
position  and  is,  in  fact,  a  strong  function  of  height  near  the 
ground.  Also  the  vapor  clouds  are  not  always  neutrally  buoyant. 
In  many  cases,  buoyant  effects  may  dominate  the  atmospheric 
transport.  For  these  reasons  more  sophisticated  models  are 
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Figure  3-2 

Temperature  Profiles  indicating  lapse  rates: 

Slightly  Unstable  (C) ,  Neutral  (Dry  Adiabatic) ,  Neutral  (D)  , 
and  Slightly  Stable  (E) . 
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Table  2.1 


Stability  Classification  by  aT/aZ 


at  in  °c/ioo  meters 

j  DEFINITION 

Less  than  -1.9 

Extremely  unstable 

-1.7  to  -1.9 

Unstable 

-1.5  to  -1.7 

Slightly  unstable 

-0.5  to  -1.5 

Neutral 

-0.5  to  1.5 

Slightly  stable 

1.5  to  4.0 

Stable 

Greater  than  4.0 

Extremely  stable 

To  accomplish  this,  2  two-dimensional  hydrodynamic  models  are 
presented.  These  allow  a  velocity  profile  to  be  present  in  the 
wind  field  and  account  for  buoyancy  terms  due  to  density 
variations . 

2.2.1  Axi-symmetric  Dispersion  Model 

In  a  calm  (  or  near  calm)  environment,  axi-symmetric  flow 
should  be  assumed  in  analyzing  the  air  entrainment  and 
dispersion  effects  of  the  well-defined  round  vapor  cloud.  In 
a  cylindrical  coordinate  system,  the  governing  equations  of  the 
time  averaged  mean  flow  variables  are: 


Conservation  of  Mass 


(9)  r  ac(eV  + 


Conservation  of  Momentum 


K b  («.  J?) 


where 

Vv  *  radial  velocity 
V2  3  axial  velocity 
p  =  pressure 

Kjt) >  Kc,  Kt  3  monentum,  concentration,  and  thermal 
dif f usivities 
p  3  density 
Po  3  ambient  density 
g  3  gravitational  constant 
T  3  temperature 
T0  3  ambient  temperature 

T  =*  T  -  Tq 

r  3  adiabatic  lapse  rate 
c  3  concentration 
t  3  time 


In  order  to  eliminate  pressure  terms  in  the  momentum 
equation,  the  vorticity  determination  is  introduced 
into  the  problem  as: 
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By  cross  differentiating  and  subtracting  the  momentum  equa¬ 
tions  (  10  )  and  (11  ),  and  using  the  conservation  of 

nass  equation,  (  9  )#  one  is  able  to  derive  an  equation 

governing  the  vorticity  as: 
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Now  two  stream  function  ^  can  be  introduced  to  satisfy 
the  continuity  equation  (  9  ): 


(16) 


v  «  I  14 

r  r  jz 


(17)  v2 


1  Jhfc 

r  ar 


Through  the  definition  of  vorticity,  the  stream  function 
will  satisfy  the  following  equation: 


(18)  3_1*  +  ihi  _  A  14  _ 
3  z2  ar2  r  ar 


This  concludes  the  system  of  equations  governing  the 
dispersion  of  an  axi-symnetric  cloud  in  quiescent 
environment. 


Quasi-Three  Dimensional  Dispersion  Model 

This  nodel  follows  a  specific  slice  of  air  as  it  moves 
downwind  under  the  influence  of  a  mean  wind  flow.  The 
hydrodynamic  equations  expressing  the  conservation  of  mass, 
momentum,  and  energy  together  with  a  Boussinesq 
approximation  for  the  atmosphere  are  written  in  rectangular 
coordinates.  In  a  coordinate  system  such  that  x  is  the 
downwind  direction,  y  is  the  crosswind  direction,  and  z  is 
the  vertical  direction,  these  equations  are: 

Conservation  of  Mass 
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where 


4  +  u-i 

at  ax 


+  v-4  +  w— 2. 


u,v,w  =  velocity  components  in  the  x,y,z 
directions,  respectively. 

If  the  assumption  is  made  that  the  downwind  velocity  is 
independent  of  position  and  unmodified  by  the  addition  of 
the  plume,  the  momentum  equation  in  the  x  direction  can  be 
omitted  (it  is  assumed  to  be  satisfied  by  the  mean 
atmosphere  motion  and  u  is  independent  of  x).  The  equations 
reduce  to: 
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By  cross-differencing  and  subtracting  the  momentum 
equations  and  using  the  conservation  of  mass  equation,  one 
is  able  to  reduce  the  number  of  equations  to  be  solved  by 
one.  Defining  a  vorticity  by: 
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and  using  the  momentum  equations  which  have  been  reduced  to 
eliminate  the  pressure  terms  yields  a  vorticity  equation: 

(31)  uiiai  +  llxal  +  Hgiil  =  2_ie.  +  7  .  rK  7n1 
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Since  y«v  =  0,  there  exists  a  strem  function  ^  such  that: 


w  =  li 
3Y 


combining  the  stream  function  definition  with  the  vorticity 
definition  resulted  in  a  Poisson  equation  relating  the  vor¬ 
ticity  to  the  stream  function. 


(33)  v2*  =  n 
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The  density  variations  are  written  in  terms  of  a 
temperature  and  concentration  variation  in  the  Boussinesq 
approximation  resulting  in  a  final  set  of  equations  as 
follows: 
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The  specific  relationship  betv/een  density,  temperature,  and 
concentration  will  depend  upon  the  equation  of  state  used. 
Typically  a  mixture  equation  is  used  where: 

(37)  p  =  (1-C)  Pq  +  Cpc 

Po  =  density  of  air 

Pc  =  density  of  added  constituent 

The  individual  densities  are  assumed  to  obey  a  perfect  gas 


2 . 3  Results 


2.3.1  Literature  Review 

A  review  of  the  literature  associated  with  the  use  of  aircraft 
fuels  and  methods  of  cleaning  up  spills  was  completed  as  a 
portion  of  this  project.  The  information  collected  during  this 
literature  survey  has  resulted  in  information  on  the  size  of 
liquid  pools  as  a  function  of  spill  size,  some  limited 
information  on  the  frequency  of  fuel  spills  in  the  Air  Force  and 
a  body  of  literature  associated  with  the  dumping  of  JP4  in  the 
atmosphere  and  the  resulting  vaporization  of  the  fuel.  This  data 
has  provided  some  insight  into  the  important  parameters  for  the 
modeling  effort.  Appendix  E  summarizes  the  specific  data 
utilized  in  the  modeling  effort  undertaken  in  this  study. 

Ciccone  and  Graves  (1)  completed  a  study  of  the  frequency  of 
spills  and  looked  at  means  of  mitigation  associated  with  these 
spills.  Table  2.2  presents  the  results  of  their  limited  study  on 
fuel  spills  at  Air  Force  facilities.  This  data  indicates  that 
spills  of  the  smallest  size  do  occur  approximately  10  to  15  times 
more  frequently  that  the  larger  spills.  However,  the  data 
indicate  that  on  the  average  one  large  spill  occurs  each  month  at 
each  Air  Force  operating  facility.  This  frequency  is  such  that 
consideration  should  be  given  to  the  associated  flammable  cloud 
problem  and  the  procedures  that  are  in  effect  in  the  vicinity  of 
the  spill.  A  knowledge  of  the  hazard  posed  by  these  spills  is 
important  in  determining  what  maintenance  procedures  or  other 
activities  can  occur  in  the  vicinity  of  an  aircraft  being  re¬ 
fueled.  The  ability  to  conduct  simultaneous  activities  with 
aircraft  being  refueled  could  significantly  decrease  the  "turn 
around"  time  for  an  aircraft. 


TABLE  2.2 


SUMMARY  OF  AVERAGE  SPILL  FREQUENCY 
(spills/month) 


No.  of 

Sample  Bases 

Small 

(0-4  gallons) 

Medium 

(5-^42  gallons) 

Large 

( > 4 2  gallons) 

23 

15.95 

4.93 

.91 

37 

13.79 

5.11 

1.22 

2.3.2  Liquid  Spill  Model  Results 
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A  series  of  calculations  were  made  for  different  size  pools  and 
different  environmental  conditions.  The  results  of  the  model 
runs  indicate  that  the  vapor  generation  rate  from  the  pool 
quickly  drops  to  a  steady  state  value  that  is  a  function  of  the 
mass  transfer  from  the  pool  and  the  heat  transfer  to  the  pool 
from  the  surroundings.  Figure  2.3  is  a  plot  of  the  calculated 
vapor  source  rate  for  a  40  gallon  JP4  spill  with  the  fuel  and 
ground  temperature  at  20  C  and  the  wind  speed  of  5  miles  per 
hour.  As  one  can  see  the  vaporization  rate  is  essentially 
constant.  This  is  caused  because  the  steady  state  solution  is  a 
balance  between  heat  into  the  pool  from  the  ground  and  heat  loss 
because  of  the  material  be  evaporated.  Steady  state  is  reached 
in  a  short  time  compared  with  the  total  time  for  the  liquid  pool 
to  be  evaporated. 

Table  2.3  shows  the  steady  state  rate  of  vapor  generation  for 
other  sizes  and  environmental  conditions.  Because  of  the  model 
assumption  of  a  constant  pool  thickness,  the  pool  diameter  is 
directly  related  to  the  volume  of  material  spilled.  For  the  40 
gallon  spills  the  pool  diameter  is  19.6  meters  and  for  the  4 
gallon  spills  the  pool  diameter  is  6.2  meters.  The  vapor  source 
rate  behaves  as  on  would  expect.  As  wind  speed  and  fuel 
temperature  are  increased  the  vapor  rate  is  increased.  As  the 
free  stream  vapor  pressure  is  increased  the  vapor  rate  is 
decreased.  Since  the  vapor  rate  is  almost  constant  the  time  to 
vaporize  the  spill  increases  as  the  vaporization  rate  decreases. 

The  results  of  the  liquid  spill  modeling  calculations  indicate 
that  the  vapor  source  rate  for  a  40  gallon  spill  can  vary  by 
almost  an  order  of  magnitude,  depending  on  the  assumed 
environmental  conditions  existing  at  the  time  of  the  spill.  This 
type  of  result  is  typical  of  releases  to  the  atmosphere  where  the 
ambient  environment  significantly  influences  the  results.  The 
other  significant  result  from  these  calculations  is  that  the 
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Spill  Volume  = 
Wind  Speed  = 
Air  Temperature  = 
Fuel  Temperature  = 
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TABLE  2.3 


SUMMARY  OF  LIQUID  SPILL  RESULTS 


SPILL 

SIZE 

(gal) 

WIND 

SPEED 

(mph) 

AIR 

TEMP. 

(°C) 

FUEL 

TEMP. 

(°C) 

Poo 

P 

VAPOR 

RATE 

(kg/sec) 

VAPOR 

TIME 

(sec) 

40 

10 

20 

20 

.755 

.4045 

285 

40 

5 

20 

20 

.755 

.2375 

450 

4 

5 

20 

20 

.755 

.0269 

425 

40 

5 

40 

20 

.755 

.3802 

300 

40 

5 

20 

0 

.755 

.1477 

770 

40 

1 

40 

20 

0 

.4145 

280 

40 

1 

20 

20 

0 

.2653 

430 

40 

5 

20 

20 

0 

.8711 

130 

por  source  rate  is  essentially  constant  and  could  be  assumed 
nstant  for  the  purpose  of  hazard  estimation. 

3.3  Spray  Model  Results 

e  question  of  droplet  distribution  resulting  from  a  break  in  a 
essurized  fuel  line  is  not  well  known.  The  size  of  droplets 


TABLE  2.4 

SUMMARY  OF  DROPLET  MODEL  RESULTS 


DROPLET 

DIAMETER 

(microns) 

RELEASE 

ANGLE 

RELEASE 
VELOCITY 
(m/s ec) 

HORIZONTAL 

DISTANCE 

(m) 

MAXIMUM 

HEIGHT 

(m) 

% 

LOSS 

200 

45 

30 

.72 

.3X5 

4.5 

1000 

45 

30 

5.25 

2.70 

5.0 

2000 

45 

30 

9.40 

4.61 

3.8 

3000 

45 

30 

12.8 

6.05 

3.1 

3000 

55 

30 

12.14 

6.78 

3.2 

3000 

35 

30 

13.1 

5.04 

2.9 

4000 

45 

30 

16.0 

7.20 

2.7 

developed  for  calculation  of  evaporation  of  JP4  jettisoned  by 
aircraft.  Because  the  problems  are  significantly  different  in 
formulation  the  only  comparison  of  significance  is  the  amount  of 
vaporization  that  occurs  in  the  first  few  meters.  por  a  1000 
micron  diameter  droplet  the  multi-component  model  predicts  that 
11%  would  be  vaporized  in  3  meters  and  13%  in  five  meters.  This 
result  is  2  to  3  times  as  large  as  that  predicted  by  the  single 
component  model  and  is  because  the  lighter  components  evaporate 
more  rapidly  than  the  model  of  this  study  predicts.  This  area 
will  require  more  investigation  during  the  proposed  phase  II 
program . 

2.3.4  Dispersion  Model  Results 

Table  2.5  shows  the  results  of  three  calculations  made  with  the 
Disco  Model  using  the  source  rate  for  40  gallon  spills  of  JP-4 . 
A  range  of  parameters  were  selected  so  that  information  would  be 
gained  on  the  range  of  expected  regions  where  the  flammable 
hazard  would  be  expected.  These  results  are  extended  to  the 
point  where  the  mean  concentration  predicted  by  the  model  is 
below  the  lower  flammable  limit  cf  1.3%  JP-4  in  air.  The  sources 
selected  for  these  calculations  represent  conditions  that  would 
cause  the  maximum  downwind  distance  to  occur.  A  stable  catagory 
of  F  was  selected  for  each  of  the  calculations  and  the  wind 
speeds  were  selected  toward  the  lower  end  of  wind  speeds. 

The  results  show  that  the  maximum  hazard  distances  occur  at  low 
wind  speeds  and  stable  conditions.  The  distances  predicted  by 
the  model  are  very  short  with  the  maximum  being  only  slightly 
over  7  meters  downwind  of  the  source.  The  vapor  rates  predicted 
by  the  droplet  or  spray  model  also  yield  similarly  short 
distances  . 

If  calculations  had  been  made  over  a  range  of  parameters  for 
stability  and  wind  speed,  the  resulting  distances  would  have  been 
significantly  shorter  than  those  presented.  There  are  however. 


TABLE  2.5 

SUMMARY  OF  DISPERSION  MODEL  RESULTS 


SOURCE 

WIND 

STABILITY 

RATE 

SPEED 

CLASS 

(KG/SEC) 

(M/SEC) 

.4145 

1.0 

F 

.8711 

5.0 

F 

.3802 

5.0 

F 

DOWNWIND 

DISTANCE 

(M) 


7.05 

1.17 

0.45 


several  unknowns  which  need  to  be  answered  to  determine  if  these 
distances  are  realistic.  The  model  in  its  present  form  assumes 
that  the  turbulent  diffusivity  of  the  atmosphere  is  constant  near 
the  ground.  There  exist  theoritical  arguments  for  this 
diffusivity  that  would  have  it  significantly  decrease  near  the 
ground.  This  decrease  in  diffusivity  could  result  in  somewhat 
longer  distances.  The  proposed  Phase  II  measurement  program 
would  address  this  question. 

In  summary  the  model  results  indicate  that  the  flammable  clould 
will  exist  only  very  near  the  source  and  that  the  possibility  of 
a  fire  being  started  more  than  a  few  meters  away  from  a  liquid 
pool  or  a  region  where  liquid  spray  exists  is  very  remote. 


3.0  CONCLUSIONS  AND  RECOMMENDATIONS 


The  results  of  this  phase  I  study  indicate  that  the  current 
procedures  and  restrictions  on  activities  during  the  refueling  of 
aircraft  may  in  many  ways  be  too  restrictive.  Models  were 
developed  for  consideration  of  the  typical  kinds  of  JP-4  releases 
that  can  occur  and  estimates  of  the  extent  of  the  flammable 
hazard  area  associated  with  these  releases  were  made.  These 
estimates  indicate  that  the  area  of  concern  is  limited  to  the 
near  proximity  of  liquid  spills  or  the  region  where  droplet 
sprays  may  reach.  The  vaporization  rate  of  JP-4  is  such  that  the 
evolved  vapor  cloud  will  be  diluted  below  the  lower  flammable 
limit  within  a  few  meters  of  the  source. 

Estimates  of  the  liquid  pool  size  indicate  that  for  a  40  gallon 
spill  the  pool  size  may  be  as  large  as  20  meters  in  diameter  so 
that  the  hazard  distance  slightly  greater  than  30  feet  from  the 
pool  center  if  a  short  dispersion  distance  is  added  for  the  vapor 
cloud.  Similarly,  the  calculations  of  spray  distance  indicate 
that  travel  distance  for  the  droplets  would  be  between  10  and  15 
meters.  This  results  is  a  hazard  distance  of  from  30  to  45  feet. 
Both  of  these  projected  distances  are  less  than  the  50  foot 
distance  presently  utilized  for  safety  purposes. 

Several  parameters  were  found  which  are  not  well  known  for  the 
releases  under  consideration.  Data  does  not  exist  that  documents 
these  parameters  in  detail  and  further  research  is  needed  to 
determine  these  parameters.  All  of  these  parameters  do  influence 
the  model  predictions  and  can  influence  the  conclusions.  The 
parameters  considered  to  be  of  significant  importance  are  as 
follows : 

•  Dispersion  parameters  (eddy  dif f usivities)  within  one  or 
two  pool  diameters  as  a  function  of  measured  meteor- 
oligical  conditions 


% 


( 


•  Droplet  size  distribution  associated  with  the  rupture  of 
refueling  lines  or  pressurized  tanks 

f  Drag  coefficient  associated  with  an  evaporating  droplet 

•  Liquid  pool  thickness  and  spread  rate 

•  Vapor  pressure  variation  over  a  liquid  pool  or  surround¬ 
ing  a  liquid  drop  moving  through  the  atmosphere 

•  Determination  of  evaporation  rate  to  determine  is  single 
component  model  is  satisfactory 

These  six  parameters  are  those  that  are  the  most  questionable  and 
could  influence  the  model  predictions  by  a  change  in  the 
parameter  value.  A  through  knowledge  of  these  parameters  will 
provide  the  assurance  that  the  model  prediction  are  reasonable. 
For  those  portions  of  the  model  were  data  was  available, 
comparisons  have  been  made  and  agreement  was  excellent.  In  other 
areas,  no  data  is  available  and  assumptions  have  been  required. 
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During  phase  I  of  this  study,  Tracer  Technologies  met  all 
technical  objectives  by  developing  two  source  rate  and  two 
atmospheric  dispersion  models  specifically  designed  to  predict 
the  hazards  associated  with  JP4  fuel  spills.  Additionally, 
limited  exercise  of  the  models,  which  indicate  hazard  distances 
much  smaller  than  those  implied  by  existing  safety  standards, 
resulted  in  the  identification  of  key  parameters  which  need 
further  definition  to  provide  models  which  have  significant 
credibility.  Phase  II  of  this  study  is  designed  to  assure  that 
all  the  key  parameters  have  been  identified  and  that  values  are 
determined  through  an  experimental  program  that  assure  these 
parameters  are  accurately  represented.  In  addition,  the 
experimental  program  would  provide  measurements  on  actual  spills 
that  would  provide  the  necessary  verification  of  the  model 
results.  No  spills  of  JP-4  have  been  conducted  to  determine  the 
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hazard  distance  and  this  necessary  effort  will  need  to  be 
undertaken  to  show  that  the  model  results  are  in  agreement  with 
actual  measurements. 

The  proposed  Phase  II  study  will  consist  of  (1)  additional  model 
calculations  to  define  the  experimental  program,  (2)  design  and 
conduct  of  JP-4  spills  at  Kirtland  Air  Force  Base,  (3) 
modification  of  the  existing  models  utilizing  the  results  from 
the  measurement  program,  (4)  merge  the  models  into  a  complete 
package,  (5)  develop  simplified  models  (nomographs  or  calculator 
models)  for  use  in  the  field  and  (6)  conduct  some  final  "hands- 
off"  validation  experiments  to  show  that  the  models  work 
properly.  This  program  will  provide  the  Air  Force  with  a  working 
documented  package  for  use  in  evaluation  of  fuel  spill  hazards. 
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Pool  is  a  small  program  that  solves  the  heat  and  mass  tansfer 
equations  for  a  well  mixed  liquid  pool  to  determine  the  vapor 
rate  as  a  function  of  time.  The  primary  dependent  variable  is  the 
pool  temperature  which  is  advanced  in  time  by  an  explicit 
differencing  scheme  utilizing  equation  as  the  basis  for  the 
differencing  scheme.  The  only  temperature  dependent  parameter  in 
the  calculation. 

The  program  consists  of  three  routines:  DRIVER,  POOL,  VPRESS. 
DRIVER  is  an  interactive  routing  that  asks  for  the  problem 
specific  input  data,  POOL  solves  the  differential  equation  and 
computes  the  time  varying  parameters  of  the  problem  and  VPRESS 
provides  the  vapor  pressure  of  JP4  as  a  function  of  temperature. 
Figure  A-l  is  a  flowchart  of  the  POOL  source  model.  The 
subroutines  of  the  POOL  model  are  relatively  straight  forward. 
Each  routine  is  briefly  described  below. 

Routine  DRIVER 

This  is  the  main  routine  of  the  model.  It  is  interactive  and 
asks  the  user  for  the  required  input  compute  the  vapor  source 
rate  for  a  specific  problem.  The  necessary  input  is  (1)  spill 
volume  in  gallons,  (2)  windspeed  in  mph  (3)  the  ambient  air 
temperature  in  °C  and  (4)  the  spilled  fuel  temperature  in  °C. 

Subroutine  POOL 

Subroutine  POOL  does  an  energy  calculation  of  the  liquid  pool  and 
computes  the  rate  at  which  vapor  is  evolved  from  the  liquid  pool. 
The  outputs  from  this  model  are  (1)  the  time  in  sec  after  the 
spill,  (2)  the  liquid  temperature,  (3)  the  evaporation  rate  in 
kg/sec  and  (4)  the  remaining  pool  mass  in  kg. 


DRIVER 


POOL 


MAIN  VARIABLES  OF  POOL  MODEL 


> 

AP  =  area  of  pool 

C  =  constant  in  evaporation  equation 
CPL  =  liquid  specific  heat 
CPV  =  vapor  specific  heat 
COND  =  heat  conduction  from  ground 
DELT  =  time  step 
DELH  =  heat  of  vaporization 

[ 

PRAD  =  pool  radius 
PMASS  =  liquid  pool  mass 
PVAP  =  vapor  pressure  of  JP4 
RAD  =  radiation  heat  transfer 

RED  »  ratio  of  ambient  vapor  pressure  to  saturation  value 
RHOV  =  vapor  densisty 
RMASS  =  mass  evaporation  rate 
SPILV  =  volume  of  spill 
SOLAR  *  solar  radiation  constant 
TAMB  =>  ambient  temperature 
TL  =*  liquid  temperature 
VAPOR  »  heat  transfer  from  vapor 
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DROP  is  a  small  interactive  program  that  solves  the  dynamic 
behavior  of  a  single  droplet  as  well  as  solving  the  intergral 
heat  and  mass  transfer  equations  around  the  drop.  Thus  the  model 
keeps  track  of  the  droplet  trajectory  as  well  as  computes  the 
changing  size  of  the  droplet  due  to  mass  transfer  from  the 
droplet.  The  primary  dependent  variables  in  the  equation  are  the 
x,y  positions  of  the  droplet,  the  corresponding  velocities  of  the 
droplet  and  the  temperature  of  the  droplet. 

The  equations  of  motion  and  the  coupled  energy  equation  are 
integrated  numerically  in  an  explicit  method  which  uses  current 
values  of  temperature,  velocity  and  position  to  compute  the  rate 
at  which  these  parameters  change  and  their  updates  each  of  the 
parameters.  The  equations  are  coupled  because  the  radius  of  the 
droplet  is  decreasing  due  to  mass  loss  and  the  drag  force  on  the 
droplet  is  a  function  of  the  radius.  In  addition  the  velocity  of 
the  droplet  affects  the  heat  and  mass  transfer  rates. 

The  program  consists  of  four  routines:  DROP,  QUAD,  EVAP,  & 
VPRESS.  DROP  is  the  main  routine  it  asks  for  input  and  computes 
the  marching  in  time.  EVAP  computes  the  mass  transfer  rate  from 
the  droplet.  QUAD  is  a  quadratic  equation  solver  and  VPRESS  is  a 
fitting  routine  that  provides  the  vapor  pressure  as  a  function  of 
temperature.  Figure  B-l  is  a  flowchart  of  the  DROP  source  model. 

The  subroutine  of  the  DROP  model  are  relatively  straight  forward. 
Each  routine  is  briefly  described  below. 

Routine  DROP 


This  is  the  main  routine  of  the  model.  It  is  interactive  and 
requires  input  from  the  user.  The  input  required  is  as  follows: 
(1)  the  initial  drop  radius,  (2)  the  angle  the  droplet  trajectory 
makes  with  the  ground,  (3)  the  initial  velocity  of  the  droplet, 
(4)  the  ambient  air  temperature  and  (4)  the  initial  temperature 
of  the  fuel  drop.  Utilizing  this  input  data  the  equations  of 


B-2 


A-  .*-k  ’  a  »  >  " 


Figure  B-l 


Flowchart  DROPLET  Source  Model 


motion  are  then  integrated  to  determine  the  dynamic  behavior  of 
the  drop. 


The  output  information  provided  by  the  model  is  (1)  the  X,Y 
position  where  X  is  the  horizontal  position  and  Y  is  a  vertical 
position  of  the  drop,  (2)  the  time  after  release  of  the  drop,  (3) 
the  U,V  velocities  corresponding  to  the  X,Y  directions,  (4)  the 
mass  of  the  droplet  and  (5)  the  evaporation  rate  from  the 
droplet . 

Subroutine  QUAD 

QUAD  is  a  small  routine  that  solves  a  second  order  algebraic 
equation  for  the  velocities  of  the  drop  at  each  time  step. 

Subroutine  EVAP 

EVAP  is  a  routine  that  integrates  the  energy  equation  for  the 
droplet  and  determines  the  evaporation  rate  for  the  droplet. 

Subroutine  VPRESS 

VPRESS  is  a  curve-fit  routine  which  returns  the  vapor  pressure  of 
JP-4  as  a  function  of  temperature. 


MAIN  VARIABLES  OF  POOL  MODEL 


ANG  *  initial  droplet  elevation  angle  of  trajectory 
ATEMP  =  ambient  temperature 

C  =  constant  in  evaporation  equation 
CPL  =  liquid  specific  heat 
CPV  =  vapor  specific  heat 
DELT  =  time  step 
DELH  =  heat  of  vaporization 
FTEMP  *  initial  droplet  temperature 
PMASS  =  droplet  mass 
PVAP  =  vapor  pressure  of  JP4 
RAD  *  radiation  heat  transfer 

RED  =*  ratio  of  ambient  vapor  pressure  to  saturation  value 
RHOV  »  vapor  densisty 
RDOT  »  mass  evaporation  rate 
TIME  *  time  of  simulation 
TA  *  ambient  temperature 
TL  «  liquid  temperature 

U  *  horizontal  velocity 

V  *  vertical  veloctiy 
VAPOR  *  heat  transfer  from  vapor 

X  *  horizontal  position 

Y  ■  vertical  position 
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For  numerical  solution,  equation  12  may  be  written  in  the 
form: 
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(c-1) 


in 

3 1 


n/r2 


-  3.  in 

p0  3r 
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\  v* . 
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The  advection  terms  in  (  C-1  )  are  solved  using  the  Crowley 

2nd  order  scheme  (  5  )  and  a  centered  second  difference  (9‘) 
is  used  on  the  horizontal  and  vertical  diffusion  terms.  The 
non-differential  term  is  evaluated  at  the  previous  time  step 
and  the  density  gradient  is  found  from  a  standard  first 
order  cell  centered  difference. 


C-2 


The  advection  terns  in  equation  (  C-l  )  are  in  the 


form: 
(C-2 ) 
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Centered  Second  Difference  Scheme 


The  diffusion  terms  in  equation  (  17  )  are  in  the 

form: 


(C-4 ) 


Let 

(C-5  )  ♦n+1»  "  (Fj  +  l/2  -  F  j  —  1/  2 )/ AX 


where 

Fj  +  l/2  =*  At  Kj+l/2  (*j  +  l  -  ♦j  )/AX 


and  the  K's  are  interface  quantities 


AXISYM  (AXI-Symmetric  Dispersion  Model) 


AXISYM  is  a  group  of  subroutines  that  uses  the  data 
produced  by  the  various  source  models  in  HAZARD  and  performs 
a  dispersion  calculation  which  is  dominated  by  the  grav¬ 
ity  spread  of  the  evolved  cloud.  Figure  c-1  is  a  sche¬ 
matic  diagram  of  the  overall  logic  used  in  the  AXISYM  pro¬ 
gram. 


The  first  subroutine  called  is  ASTART  which  relates 
the  variables  from  HAZARD  with  the  variables  used  in 
AXISYM.  AXISYM  than  enters  a  loop  which  advances  the  simu¬ 
lation  in  time.  Inside  of  the  loop  at  each  time-step,  the 
following  routines  are  called:  ITLAPL,  VEL,  SETDIF,  FIXUP, 
DELTAT,  0UT1,  and  REZONE.  These  routines  calculate  the 
stream  function,  velocities,  diffusities,  and  then  update 
the  vorticities  in  each  zone.  A  new  time-step  is  calcu¬ 
lated  and  data  is  printed  which  documents  the  changing 
nature  of  the  cloud. 

The  cyclic  behavior  continues  until  the  cloud  con¬ 
centration  is  below  the  flammable  limit  and/or  the  lOOppm 
HjS  toxic  limit.  A  brief  description  of  each  subroutine 

is  described  in  the  following  section.  Due  to  the  simi- 
* 

larities  between  the  routines  used  for  DISCO  and  AXISYM, 
the  reader  is  referred  to  the  flow  diagram  of  the  cor¬ 
responding  routine  in  DISCO. 

Subroutines  of  AXISYM 


Subroutine  ASTART 

ASTART  initializes  the  variables  and  arrays  for  the 
beginning  of  the  calculation.  Several  constants  and  para¬ 
meters  are  also  calculated  in  this  subroutine.  (see  Fig¬ 
ure  D-2  ?  INPUT1) 


Subroutine  ITLAPL 

ITLAPL  iteratively  solves  the  Poisson  equation  for 
the  stream  function  using  an  over-relaxation  method. 

ITLAPL  calls  source  to  determine  the  boundary  conditions 
on  the  stream  function  that  will  allow  the  correct  gas 
volume  to  flow  into  the  grid. 

Subroutine  VEL 

VEL  calculates  the  horizontal  and  vertical  velocity 
throughout  the  grid  from  the  stream-function  array. 

Subroutine  SETDIF 

SETDIP  uses  the  local  velocity  and  temperature  strati¬ 
fication  in  each  zone  to  determine  the  horizontal  and  verti¬ 
cal  dif fusivities  throughout  the  grid.  These  dif fusivities 
are  updated  on  each  cycle  of  the  calculation.  (see  Figure  d-3 
SETEV) 

Subroutine  FIXUP 

FIXUP  is  the  subroutine  where  the  conservation  equations 
for  mass,  energy,  and  vorticity  are  updated.  These  updates 
are  done  by  considering  diffusion,  advection,  and  other  terms 
such  as  the  buoyancy  and  source  terms.  The  advective  updates 
are  accomplished  in  subroutine  DIFF.  For  the  concentration 
and  temperature,  a  variable  transformation  is  accomplished 
by  multiplying  by  the  local  radius  before  calling'CROWLY 
and  DIFF.  (see  Figure  D-5  ;  UPDATE) 

Subroutine  0UT1 

0UT1  is  the  output  routine  where  information  on  the 
cloud  size  and  shape  is  provided  to  the  user. 


Subroutine  REZONE 

REZONE  is  a  routine  that  provides  the  capability  to 
enlarge  the  region  of  interest.  When  the  disturbance 
approaches  the  boundary  of  the  grid,  then  the  horizontal 
or  vertical  zone  size  is  doubled  and  new  ambient  regions 
are  added  around  the  disturbance.  (see  Figure  D-8  ) 
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Main  Program  Variables  for  AXISYM 
concentration  array  (v/v) 

change  in  concentration  array  during  time  step  (v/v) 
initial  concentration  (v/v) 
concentration  of  source  gas 
average  specific  heat  of  gas  (J-KG/°K) 
time  step  (sec) 
initial  time  step  (sec) 
horizontal  zone  size  (m) 

-  array  of  vertical  zone  sizes 

-  vertical  height  of  bottom  zone  (M) 

2 

-  radial  diffusivity  array  (M  /sec) 

2 

-  vertical  diffusivity  array  (M  /sec) 

-  gravity  (9.8  M/sec) 

-  adiabatic  lapse  rate  (°K/100  meters) 

-  surface,  heat  transfer  coefficient 

-  cycle  number 

-  initialization  flag  for  ITLAPL 

-  initialization  flag  for  SETDIF 

-  IZONE  is  initialization  flag  for  REZONE 

-  total  number  of  cycles  in  problem 

-  number  of  zones 

-  NX-1 

-  NX-2 

-  NX+1 

-  number  of  vertical  zones 

-  NY-1 

-  NY-2 

-  NY+1 

-  previous  time  setp  (sec) 

-  stream  function  (M^/sec) 

-  pressure  divided  by  specific  gas  constant 

-  1/ ( 2* AX) 

-  1/DX 
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RX 

RXMIDI 

SLOPE 
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-  1/ (DX) 2 

-  array  of  radial  distance  to  inner  zone  boundary 

-  array  of  1/RXMID 

-  slope  of  temperature  profile 

-  specific  gravity  of  gas 

-  temperature  array  (°K) 

-  initial  ground  level  temperature  (°K) 

-  change  in  temperature  array 

-  adiabatic  temperature  profile  (°K) 

r  perturbation  in  temperature  profile  form  ground  (°K) 

-  time  (sec) 

-  initial  temperature  of  gas  (°K) 

-  radial  velocity  (M/sec) 

-  vertical  velocity  (M/sec) 

-  vorticity  array  (1/sec) 

-  changej  in  vorticity  array  during  time  setp  (1/sec) 

i 

-  initial  vorticity  array  (1/sec) 

-  array  of  heights  to  bottom  of  zone 

-  ratio  of  vertical  zone  growth 


Equations  (34  ) ,  (  35  ),  and  (  36  ')  are  solved  by 

nunberial  integration  with  the  integration  stepping  forward 
in  terns  of  downwind  distance  x.  These  equations  are 
updated  and  the  new  values  are  saved  for  use  during  the  next 
tine  or  calculational  cycle.  Equation  (33  )  is  then 

solved  and  new  velocities  are  calculated  from  the  gradient 
of  the  stream  function.  The  distance  interval,  &x,  is 
selected  to  satisfy  numerical  stability  criteria  associated 
with  the  advection  or  horizontal  diffusion  terns  of  each 
equation. 

Again,  as  in  the  axi-symmetric  model,  the  advection 
terns  are  solved  using  the  Crowley  2nd  Order  Scheme.  The 
horizontal  diffusion  term  is  solved  using  the  centered 
second  difference  scheme  and  the  vertical  diffusion  is 
solved  from  an  implicit  scheme  discussed  below. 


Implicit  Scheme  for  Vertical  Diffusion 

The  implicit  method  used  requires  that  certain 
intermediate  quantities  be  stored  while  looping  upwards  and 
then  the  *  quantity  is  updated  looping  downwards. 
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Since  this  method  requires  immediate  updating  of  the  * 
quantity,  this  is  solved  last.  The  other  terms  in  this 
expression  are  given  by: 
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Stream  Function  Calculation 


Equation  (  33  )  for  the  stream  function  is  solved  by 

means  of  Fourier  transforms  resulting  in  a  direct  solution 
of  the  differential  equation.  A  second  order  finite 
difference  approximation  to  the  Poisson  equation  a2\|i  =  is 
obtained  by  replacing  the  second  derivative  operator  by  a 
centered  second  difference  operator. 


(  D-9  ) 
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i  =  1/2  . 
j  =  2 1  •  • 


where 

6y  *ij  =  +i+l,j  “  2*ij  +  *i-l,j 

and 


6 2  +ij  =  *if j+1  ~  + 

Boundary  conditions  are  imposed  as  follows: 


,  I 
J-l 


At  the  bottom  of  the  mesh, 

^  ^  ^  =  <x^  r  i  s  1|  •••)  It 

At  the  top  of  the  mesh, 

♦  ®  8^  #  i  =  1  »2 ,  •  •  •  *  I 

The  cyclic  boundary  conditions  in  the  horizontal  area, 

^o,j  ^I»j  and  ^l,j  *  ^I+lrj  ,  j  «  2,  J-l. 

Ue  introduce  an  orthonormal  base  set  of  functions 
having  cyclic  properties  on  the  index  i: 
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where  I  is  even 

k  =  1  f  2  f  •  •  •  f  1  • 


These  are  the  finite  Fourier  functions,  which  have  the 
properties, 

I 

(D-10)  \  wik  wij,  =  ®kji 
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and  the  analogous  cyclic  boundary  conditions  are  valid  in 
the  horizontal.  They  also  have  the  property  that  they  are 
eigenfunctions  of  the  central  second  difference  operator 
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where  *k  =  2  sin  wk/I .  These 
functions  on  the  interval  i  = 
arbitrary  function  f^  on  this 


functions  are  complete 

1,2,  . ..,  I.  Consequently,  an 

space  can  be  represented 
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Vie  are  now  ready  to  consider  equation  (  d-9  '  )  from  the 
point  of  view  of  Fourier  transformation.  The  vorticity  and 
stream  function  are  represented  as  Fourier  series  as 
follows: 
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(1)  The  vorticity  and  the  top  and  bottom  boundary 
values  of  the  stream  function  are  subjected  to  Fourier 
transformation  to  obtain 
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and 
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(2)  The  Fourier  components  of  the  stream  function  are 
obtained  by  solving  the  tridiagonal  system  of  equations, 
equation  ^  d-15^  f°r  aji* 

(3)  The  stream  function  itself  is  obtained  by  Fourier 
synthesis 


ij  “tl1  aU  Wi t 


The  quantity  I  must  be  even.  In  order  to  take  maximum 
advantage  of  the  efficiency  of  the  Fast  Fourier  Transform, 
the  quantity  I  should  also  be  a  power  of  2. 


D-6 


we  obtain: 


In  equation  (D-15  >  the  value  of  the  wave  number,  *, 
appears  only  parametrically.  For  each  value  of  t  there  is  a 
tridiagonal  equation  having  fixed  values  at  the  end  points 
of  the  j-interval.  m 

Summarizing  the  procedure  for  obtaining  the  direct 
solution  of  the  Poisson  equation,  equation  (  d-9  ),  by 
Fourier  transform: 


Program  DISCO 

The  DISCO  model  is  a  two-dimensional  hydrodynamics  atmospheric 
dispersion  code  that  ignores  all  interactions  in  the  downwind 
direction.  The  model  solves  the  conservation  equations  for 
energy,  mass,  and  momentum  as  discussed  in  Section  2. 


This  model  calculates  the  movement  of  a  cloud  in  three 
dimensions  by  marching  a  slice  of  air  downwind  from  the  source. 

This  two-dimensional  plane  is  perpendicular  to  the  wind  direction 
and  is  advected  at  the  ambient  wind  velocity. 

There  are  a  total  of  20  subroutines  and  a  main  program 
within  the  DISCO  model.  The  organization  of  all  subroutines 
are  shown  in  Figure  D-l.  Appropriate  comment  cards  have  been 
inserted  throughout  the  code  in  order  to  provide  an  easily 
understandable  program. 

The  DISCO  main  program  is  the  driving  routine'  that  contains 
the  primary  loop  in  the  program.  The  first  subroutine  called 
is  INPUT1  which  intializes  several  variables.  The  input  parameters 
are  obtained  from  dLINK  and  are  passed  through  the  common  ZZBILL. 
Then  the  routine  SETUP  is  called.  This  is  a  subroutine  to  DISCO 
and  is  not  the  same  as  the  interactive  SETUP  program.  SETUP 
initializes  the  coefficients  to  be  used  by  the  Fast  Fourier  Transfer 
routine. 

♦ 

The  main  loop  of  the  program  is  entered  next.  First  the 
Poisson  equation  is  solved  and  the  velocities  in  the  plane  are 
calculated.  Then  the  eddy-dif fusivities  are  computed  based  on 
the  local  stability  and  velocities.  A  time  step  is  calculated 
for  the  computational  cycle  based  on  numerical  stability  criteria. 
Then  the  concentration,  temperature,  and  vorticity  arrays  are 
updated  by  transporting  the  values  with  the  local  velocities  and 
diffusing  the  values  with  the  calculated  diffusivities. 


The  output  routine  is  entered  and  the  cloud  width  to 
several  concentrations  are  printed.  The  time  (actually 
distance)  is  updated  and  the  grid  is  checked  to  make  sure 
the  concentration  of  the  pollutants  is  negligible  at  the 
boundaries.  If  the  concentration  is  significant,  the  grid 
will  be  rezoned  in  either  the  horizontal  or  vertical  direc¬ 
tion  or  both. 

The  final  routine  SAVER  checks  for  completion  of  the 
run.  Also  if  the  changes  in  the  velocities  are  small,  SAVER 
skips  the  calculation  of  velocities,  diffusities,  and  the 
time  step. 

Subroutines  of  DISCO 

Subroutine  INPUT1 

INPUT 1  initializes  all  of  the  variables  and  indices  to 
be  used  for  the  problem.  The  input  is  obtained  from  DLINK 
through  common  ZZBILL.  Also  the  grid  dimensions  and  ambient 
profiles  of  pressure  and  temperature  are  generated.  (see 
Figure  D-2  ) 

Subroutine  SETUP 

This  routine  sets  up  coefficients  for  the  Fast  Fourier 
Transfer  routines.  See  Appendix  D  for  details.  All  FFT  rou¬ 
tines  are  available  in  the  literature  and  are  not  described 
in  detail  here. 

Subroutine  LAP LAC 

This  routine  controls  the  solution  of  the  Poisson's 
equation.  See  Appendix  £>-. 

Subroutine  FFANL 

Fast  Fourier  Analyzer.  See  Appendix  D. 

Subroutine  GAUSS 

Gauss  Elimination  in  the  vertical.  See  Appendix  D. 


SUBROUTINE  INPUT  1 


Subroutine  FFSYN 

Fast  Fourier  Synthesizer.  See  Appendix  D. 

Subroutine  XYPOIS 

FFT  routine  used  to  call  FFANL,  GAUSS,  and  FFSYN. 
Subroutine  COOTUK 

Cooley-Tukey  Fast  Fourier  Transform  Routine 
Subroutine  VELOC 

VELOC  calculates  the  velocity  field  based  on  the  stream 
function  psi.  See  Appendix  D . 

Subroutine  SETEV 

SETEV  calculates  the  local  eddy-diffusivities  based  on 
atmospheric  stability  and  wind  velocities.  The  data  used  is 
empirical  from  Smith  and  Howard.  (see  Figure  D-3  ) 

Subroutine  TIMSTP 

The  routine  TIMSTP  calculates  the  internal  time  step 
based  on  numerical  stability  criteria.  Actually,  time  is 
a  distance  in  the  computer  code.  (see  Figure  d-4  ) 

Subroutine  UPDATE 

*  _  -  - 

This  routine  calculates  the  changes  in  the  scalar  vari¬ 
ables  due  to  advection,  diffusion,  and  heat  transfer  effects. 
Then  this  routine  updates  the  old  valves  of  temperature,  con¬ 
centration  and  vorticity  to  find  the  new  valves  at  this  down¬ 
wind  distance.  (see  Figure  d-5  ) 

Subroutine  0RDER2 

0RDER2  calculates  the  changes  in  the  temperature 
concentration  and  vorticity  arrays  due  to  advection.  These 
advection  terms  are  solved  utilizing  a  Crowley  second  order 
scheme.  (see  Figure  D-6  ) 


Subroutine  SETEV 


SUBROUTINE  TIMSTP 


SUBROUTINE  UPDATE 


Zero  storage 
arrays 


Calculate  the 
buoyancy  term 


Calculate  the 
surface  heating  term 


Calculate  the 
advection  terms 
CALL  ORDER  2 


Calculate  the 
diffusion  terms 
CALL  DIFFUS 


Update  the 
arrays 


Correct  numerical 
errors 


Return 


Figure  D-5 


SUBROUTINE  ORDER  2 


Subroutine  DIFFUS 


DIFFUS  calculates  the  diffusion  terms  using  the  eddy 
viscosities  computed  in  SETEV.  The  horizontal  terms  are 
calculated  using  a  standard  second  order  centered  dif¬ 
ference  scheme.  The  vertical  terms  are  solved  using  an 
implicit  method  (see  Figure  0-7  ) • 

Subroutine  OUTPUT 

OUTPUT  prints  out  the  output  of  cloud  width  at  a 
downwind  distance.  The  width  is  found  to  the  lower  flam¬ 
mable  limit,  the  stoichiometric  limit,  and  to  two  toxicity 
limits  of  if  required.  Also  the  cloud  height  and  the 
maximum  concentration  at  each  downwind  distance  are  printed 

Subroutine  RE 2 ONE 

If  s  significant  concentration  of  pollutant  nears  the 
edge  of  computational  grid,  the  grid  is  rezoned  to  contain 
the  entire  problem  within  the  grid.  REZONE  can  change  the 
grid  scale  in  either  vertical  or  horizontal  direction  or 
both.  (see  Figure  D-ff  ) 

Subroutine  SAVER 

SAVER  is  the  routine  that  determines  if  the  problem 
is  completed.  This  involves  testing  on  whether  the  maxi- 
mum  concentration  is  less  than  the  lower  flammable  limit. 

If  H 2S  is  considered  in  the  problem,  then  the  maximum  con¬ 
centration  is  also  compared  to  the  lOOppm  limit  of  H2S. 

In  addition,  for  situations  where  the  velocity  field  is  not 
changing,  the  SAVER  routine  skips  LAPLAC,  SETEV,  and  TIMSTP 
(see  Figure  D-9  ) 


Subroutine  GETMAX 

GETMAX  finds  the  maximum  valve  of  the  incoming  array. 
Subroutine  GETMIN 

GETMIN  finds  the  minumum  valve  of  the  incoming  array 


SUBROUTINE  DIFFUS 


SUBROUTINE  RE ZONE 


KTYPE 


1  or  3 


Horizontal  rezone — 
Reinitialize  grid  parameters 


Left  half  of  grid — 
average  old  valves  and 
store  in  quarter  of  grid 


Setup  new  section  of 
left  half  of  grid 


Right  half  of  grid 
average  old  valves  and 
store  in  quarter  of  grid 


Setup  new  section  of 
right  half  of  grid 


Update  C,T  &  W  arrays 


KTYPE  ^ 

\  ^2  or  3 


Vertical  rezone  -  - 
Reinitialize  grid 
parameters 


Average  old  valves  and 
store  in  bottom  half  of 


Setup  top  section  of 


Update  C,T  &  W  arrays 


Setup  PFT  coefficients 
CALL  SETUP 


c 

Cl 

CINIT 

CLIM 

CMAX 

CPMETH 

DIFDAT 

DIFNOW 

DT 

DX 

DY 

DYCUM 

EVX 

EVY 

HEATTC 

IONCE 

MTEMP 

NX 

NX2 

NXM1 

NXM2 

NXP1 

NY 

NYM1 

NYM2 

NYP1 

P 

PRMETH 

R2DX 

RDX 

RDX2 


MAIN  Program  Variables  for  DISCO  Model 

-  pollutant  array  (volume  concentration) 

-  stores  one  cycle  changes  in  pollutant  array 

-  initial  pollutant  concentration 

-  array  containing  flammable,  toxicity  and  stoich¬ 
iometric  limits 

-  maximum  pollutant  concentration  for  this  cycle 

-  heat  capacity  of  pollutant  (c/gk°) 

-  empirical  diffusion  coefficients 

-  interpolated  diffusion  coefficients  for  current  grid 

-  distance  step  for  this  cycle  (meters) 

-  grid  spacing  in  the  horizontal  direction  (meters) 

-  grid  spacing  in  the  vertical  direction-variable  (meters) 

-  cumulative  heights  (meters) 

-  horizontal  eddy  diffusivity 

-  vertical  eddy  diffusivity 

-  heat  transfer  coefficient  (J/kg°c) 

-  flag  for  SETEV  if  grid  re zoned 

-  adiabatic  lapse  rate  (°c/m) 

-  grid  zones  in  x  direction 

-  NX  divided  by  two 

-  NX  minus  one 

-  NX  minus  two 

* 

-  NX  plus  one 

-  grid  zones  in  vertical  direction 

-  NY  minus  one 

-  NY  minus  two 

-  NY  plus  one 

-  stream  function  array 

-  pressure  divided  by  gas  constant  for  pollutant 

-  one  half  divided  by  DX 

-  one  divided  by  DX 

-  one  divided  by  DX  squared 


RR 

RTEMP 

SI 

SL 

SPGRAV 

T 

TI 

TAMB 

TIME 

TKELV 

TLNG 

TSURF 

TT 

0 

01 

UFIX 

UU 

V 
W 
W 
W1 
WTO  I 

Y 

YBOTM 

YRATIO 


-  coefficients  used  in  FFT 

-  atmospheric  lapse  rate 

-  coefficient  used  in  FFT 

-  coefficient  used  in  FFT 

-  specific  gravity  of  pollutant 

-  temperature  array  (°C) 

-  stores  one  cycle  changes  in  temperature  array 

-  ambient  temperatures  (°C) 

-  distance  downwind  (meters) 

-  temperature  Kelvin  (273.16°) 

-  initial  temperature  of  pollutant  (°K) 

-  surface  temperature  (°C) 

-  coefficient  used  in  FFT 

-  horizontal  velocity  array  (m/s) 

-  stores  old  values  of  u  velocity 

-  initial  wind  velocity  array  (m/s) 

-  coefficient  used  in  FFT 

-  vertical  velocity  array  (m/s) 

-  coefficient  used  in  FFT 

-  vorticity  array 

-  stores  one  cycle  changes  in  vorticity  array 

-  gravitational  constant  multiplied  by  TKELV 

-  heights  at  center  of  each  vertical  zone  (m) 

-  height  of  bottom  vertical  zone  (m) 

-  ratio  of  heights  between  each  vertical  zone 


Calculation  of  Droplet  Size  Distribution 

An  analysis  of  the  JP4  fuel  venting  from  Aircraft  inflight  (El) 
indicates  that  the  expected  droplet  size  would  range  from  0.5 
micron  to  2100  microns.  The  larger  number  is  the  maximum  stable 
drop  size  for  a  free  falling  drop  in  the  atmosphere  whereas  the 
small  drop  is  estimated  to  be  formed  from  the  shearing  between 
the  liquid  stream  and  the  air  velocity  caused  by  the  aircraft 
forward  speed. 

The  maximum  velocity  expected  from  a  fuel  line  break  would  be 
appr oximatedly  30  m/sec  which  is  significantly  less  than  aircraft 
speed.  This  would  suggest  that  the  sprayed  droplets  would  tend 
to  be  somewhat  larger  that  would  be  expected  from  the  fuel 
dumping  from  an  inflight  aircraft.  In  actual  situations,  a 
single  droplet  size  would  not  be  produced  but  a  distribution 
could  be  expected.  In  addition,  a  significant  time  would  be 
required  for  oversized  droplets  to  become  unstable  and  "break-up" 
into  smaller  droplets.  The  fuel  dump  experiments  (El)  resulted 
in  a  measured  droplet  distribution  with  the  median  drop  diameter 
of  several  100  microns.  For  the  problem  of  spray  from  a  failed 
pressurized  line,  the  sizes  are  expected  to  be  larger  and  result 
in  droplet  sizes  of  several  thousands  microns.  for  our  nominal 
calculations  a  3000  micron  droplet  diameter  was  used. 

Calculation  of  Drag  Coefficient  for  Spherical  Droplets 

Figure  E-l  presents  data  for  the  drag  coefficient  Cp  of  a  sphere 
moving  through  air  as  a  function  of  the  Reynolds  number.  For 
this  study,  an  approximate  Reynolds  number  may  be  calculated  as: 


O')  1*3 


Coefficient  as  a  Function  of  Reynolds 
Number  from  (E2) 


L."l  ■  <  l  I  I  JH'.  I  1.^ 


i 


where 


y  =  air  viscosity 
D  =  droplet  diameter 
Pa=  air  density 


The  velocity  of  the  droplet  may  be  calculated  using  Bernoulli's 
equation  as  indicated  below: 


AP  =  1  V2 
P  2  g 


C 

For  a  fueling  operation,  the  fuel  line  pressure  is  approximately 
50  psig.  The  velocity  associated  with  a  droplet  emitted  from  a 
hole  in  the  line  is  calculated  to  be  96  ft/sec.  (  s29  m/sec). 
Figure  E-2  shows  the  drag  coefficient,  Cp  ,  as  a  function  of 
droplet  diameter  and  wind  speed  assuming  standard  atmospheric 
conditions.  Reynolds  numbers  vary  from  approximately  5  (U=  1 

^  m/sec,  D=500  microns)  to  9500  (U  =  50  m/sec,  D=2000  microns). 

For  the  purposes  of  the  phase  I  study,  the  CD  was  assumed 
constant  at  a  value  of  approximately  0.4  ie,  the  approximate 
value  for  droplet  diameters  greater  than  500  at  a  velocity  of 
30  m/sec.  During  phase  II,  calculations  will  be  performed  to 
determine  if  further  definition  of  this  variable  is  required. 
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Calculation  of  Turbulent  Eddy  Dif fusivitics 

For  the  model,  an  ambient  turbulence  model  according  to 
Hanna  (E3)  and  correlated  by  Smith  and  Howard  (E4)  was 
adopted.  Hanna  has  suggested  a  form  for  the  equation 
expressing  the  local  diffusivity  which  can  be  expressed  in 
any  of  the  following  equivalent  forms, 

(E-l)  Ka  =  0.09  0wkra  “l 
(E-2)  Ka  =■  0.06  e1/3  k^"1 
( E—  3 )  Ka  *  0.35  aw4  e“x 

where  aw  i-3  the  standard  deviation  of  the  vertical  turbu¬ 
lence,  c  is  the  rate  of  dissipation  of  turbulent  energy,  and 
k„  is  the  wave  number  where  the  vertical  turbulent  energy 
spectrum  is  a  maximum.  Hanna  evaluated  the  constants  given 
in  these  equations  through  the  use  of  observational  data. 

The  ambient  turbulence  eddy  diffusivity  Ka  was  evaluated 
alternately  by  equations  with  arbitrary  constants  and  an 
equation  relating  the  atmospheric  momentum  flux  ,to  the  mean 
velocity  shear.  The  relation  between  the  two  yielded  the 
values  of  the  constants  given. 

Numerous  measurements  of  aw  have  also  been  made  under 
a  variety  of  conditions.  Alternately,  Ka  was  related  to  a 
turbulence  scale  length  (t)  by  Taylor,  Warner  and 
Bacon  (es)*s 


Because ‘of  the  ready  availability  of  observational 
data,  this  study  uses  a  solution  technique  where  Ka  is 
calculated  as 

(E-5)  Ka  *  0.45  <jw t  *  0.45  u<je£ 

since  „w  is  equal  to  u<je  where  u  is  the  local  velocity 
and  o  is  equivalent  to  a  measurement  of  wind  vane 
fluctuations  in  radians.  The  values  of  0z  and  t  are 
functions  of  height  and  atmospheric  stabilities  and  the 
values  used  in  this  study  are  summarized  in  Tables  E-l  and 
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TABLE  E-l 


Variations  in 


TABLE  E-2 


Turbulence  Scale  Lengths  (t) 


i(m) 


BCD 


P  G 


18  15  12  10  8  7  6 

30  25  21  18  16  14  12 

41  34  29  25  22  20  17 

62  52  44  39  35  31  27 

84  71  60  52  48  43  37 

105  85  74  64  60  54  46 


diffusivity  is  relatively  enhanced  over  the  vertical,  while 
in  an  unstable  atmosphere  the  opposite  is  the  case,  although 
to  a  lesser  degree.  In  Table  e-3  we  give  the  ratio  of 
horizontal  to  vertical  diffusivity  as  a  function  of  Pasquill 
category  stability  classes  according  to  Lantz  (E6  ). 

Based  on  field  data  gathered  over  four  years  by  SAI 
(  E7) ,  values  of  KafX/KafZ  of  1.0,  10.0,  and  25.0  have 
been  used  in  the  calculations  presented  for  stability 
classes  D,  E,  and  F.  The  values  are  conservative  compared 
to  those  of  Lantz. 

TABLE  E-3 

Ratio  of  Horizontal  and  Vertical 
Dif fusivities  Versus  Stability  Class 


Stability 


Ratio 


SOURCE  RATE  SPECIFICATION  FOR  DISPERSION  MODELS 

The  axi-symmetric  dispersion  model  requires  a  source  rate  be 
specified  that  is  both  time  and  position  dependent.  For  the 
pipe  rupture  the  source  codes  provide  the  time  dependent 
release  rate  and  the  release  is  positioned  at 


the  break  in  the  first  grid  cell.  The  tank  failure  pool 
spread  models  give  both  the  pool  position  and  vaporization 
rate  as  functions  of  time  for  direct  input  into  the 
dispersion  model. 

The  quasi-three  dimensional  dispersion  model  requires  a 
steady  state  source  rate  be  specified.  The  source  rate  data 
provided  by  the  source  modesl  is,  therefore,  converted  into 
an  equivalent  steady  source  using  an  averaging  scheme 
consistent  with  the  fundamental  assumptions  involved  in  the 
quasi-three  dimensional  dispersion  model.  The  grid  is 
positioned  at  the  break  for  pipeline  failures  and  at  the 
center  of  the  pool  in  the  liquid  spills. 


Source  Model  Parameters 


Ciccone  and  Graves  (E8),  in  a  study  of  fuel  tank  leak  classifica¬ 
tion,  measured  the  surface  area  of  spill  as  a  function  of  the 
spill  size.  An  analysis  was  made  of  the  data  presented  in  the 
report  and  it  was  concluded  that  the  thickness  of  the  liquid  pool 
would  be  approximately  0.0005  meters  (.5mm).  This  number  was 
used  to  estimate  the  final  pool  spill  size  in  the  model. 

The  other  information  obtained  from  the  literature  was  the  vapor¬ 
ization  rate  for  a  droplet.  Mackay  and  Matsugu  (E9)  provide  the 
basic  theory  for  liquid  spills  of  hydrocarbons  on  land.  However, 
one  of  the  unknowns  in  that  theory  is  the  average  vapor  partial 
pressure  above  the  liquid  pool.  The  maximum  evaporation  rate  can 
be  determined  by  assuming  that  the  vapor  partial  pressure  in  the 
atmosphere  is  zero.  Dawbarn,  Nutt  and  Pender  (E10)  have  made 
measurements  on  the  vaporization  rate  of  small  droplets  in  a  wind 
tunnel.  This  information  on  the  droplet  vaporization  rate  was 
compared  with  the  theory  and  it  was  determined  that  a  nominal 
partial  pressure  of  .755  times  the  vapor  pressure  provided  the 
best  fit  to  the  data.  Figure  E-3  provides  a  comparison  of  the 
measured  data  and  the  fit.  This  fit  or  P  =0  was  utilized  for 

OO 

calculations  presented  in  this  report  giving  a  range  of  hazard 
distances  as  a  funtion  of  this  parameter. 
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Figure  E-3 

Droplet  Evaporation  Rate  for  a  500  Micron  Radius 
Droplet  as  a  Function  of  Time  (U=1Q  mph,  T=20°C) 
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Spill  Volume  (cu  m)  =  .1514 

Wind  Speed  (m/sec)  =  2.23 

Air  Temperature  (c)  =  20.00 

Fuel  Temperature  (c)  =  20.00 

Pool  Radius  (m)  =  9.82 


TIME 

SEC 

.  0000E+0 1 
.  5000E+01 
.  1000E+02 
.  1 500E+02 
.  2000E+02 
.  2500E+02 
.  3000E+02 
.  3 500E+02 
.  4000E+02 
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.  2250E+03 


TEMP 

C 

.  2932E+03 
.2921 E+03 
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VAPOR  RATE 
KG/SEC 
.  0000E+0 1 
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POOL  MASS 
KG 
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DROPLET  MODEL 
DROP  RADIUS  =  .  1500E+04 

RELEASE  ANGLE  =  45.000 

RELEASE  VELOCITY  =  .3000E+02 

AIR  TEMPERATURE  =  20.0 
FUEL  TEMPERATURE  =  20.0 


.  9044E  +  00 
.  1693E+01 
.  2392E+01 
.  3019E+01 
.  3589E+0 1 
.  4 1 10E+0 1 
.  4590E+01 
.  5036E+0 1 
.  5451E+01 
.  584 1E+01 
.  6207E+0 1 
.  6552E+01 
.  6879E+01 
.  7190E+01 
.  7485E+01 
.  776SE+01 
.8037E+01 
.  8296E+01 
. 854  3E+0 1 
. 8782E+01 
. 901 IE+01 
. 9232E+01 
. 9446E+0 1 
. 9652E+01 
.9851 E+0 1 
. 1004E+02 
. 1023E+02 
. 1041E+02 
. 1059E+02 
. 1076E+02 
. 1093E+02 
. 1109E+02 
. 1125E+02 
. 1140E+02 
. 1155E+02 
.  U70E+02 
. 1185E+02 
.  U99E+02 
. 1212E+02 
.  1226E+02 
. 1239E+02 
.  1252E+02 
. 1264E+02 
.  1277E+02 
.  1289E+02 


.  8833E+00 
.  1634E+01 
. 2281E+01 
. 2844E+01 
.  3338E+0 1 
. 3772E+01 
.  4 1 54E+01 
. 4491E+01 
. 4787E+01 
.  5046E+01 
. 5272E+01 
.  5465E+0 1 
.  5629E+01 
.  5765E+01 
. 587  3E+0 1 
. 5956E+0 1 
. 6014E+01 
. 6046E+01 
. 6054E+0 1 
. 6038E+01 
. 5997E+01 
. 5932E+01 
. 5844E+01 
. 5732E+01 
. 5599E+01 
. 5445E+01 
. 5270E+01 
. 5077E+01 
. 4865E+01 
. 4637E+01 
. 4394E+01 
. 4136E+01 
.  3865E+0 1 
.  3583E+01 
. 3289E+01 
.  2985E+01 
. 2673E+01 
. 2352E+0I 
.  2024E+01 
.  1 689E+01 
.  1349E+01 
.  1003E+01 
. 6524E+00 
.  2978E+00 
- . 6056E-01 


TIME 

.  5000E-0 1 
.  1000E+00 
.  1500E+00 
.  2000E+00 
.  2500E+00 
.  3000E+00 
.  3500E+00 
.  4000E+00 
.  4500E+00 
.  5000E+00 
.  5500E+00 
.  6000E+00 
.  6500E+00 
.  7000E+00 
.  7  500E+00 
.  8000E+00 
.  8  500E+00 
.  9000E+00 
. 9500E+00 
. 1000E+01 
. 1050E+01 
. 1100E+01 
. 1 1 50E+0 1 
. 1200E+01 
. 1250E+01 
. 1 300E+0 1 
. 1350E+01 
. 1400E+01 
. 1450E+01 
. 1500E+01 
. 1550E+01 
. 1600E+01 
. 1650E+01 
. 1700E+01 
.  1 7  50E+0 1 
. 1800E+01 
. 1850E+01 
. 1900E+01 
.  1950E+01 
.  2000E+01 
.  2050E+0 1 
. 2100E+01 
. 2 150E+01 
.  2200E+01 
. 2250E+01 


. 2870E+01 
.  2804E+01 
.  2741E+01 
.  268  IE+01 
.  2624E+01 
. 2569E+01 
. 2516E+01 
.  2465E+01 
.  241 6E+01 


- . 6253E+0 1 
- . 641 5E+01 
- .  6561E+01 
- .  6692E+01 
-  .6810E+01 
- . 6915E+01 
- .  7009E+01 
- . 7092E+01 
- . 7166E+01 


. 1046E-04 
. 1045E-04 
. 1044E-04 
. 1044E-04 
. 1043E-04 
.  1042E-04 
.  1041E-04 
.  1041E-04 
.  1040E-04 
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u 
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MASS 

RDOT 

. 1809E+02 

. 1767E+02 

. 1074E-04 

.  3996E-06 

. 1577E+02 

. 1501E+02 

. 1072E-04 

. 3525E-06 

. 1398E+02 

.  1294E+02 

.  1071E-04 

.  3 1 53E-06 

. 1255E+02 

. 1126E+02 

. 1 069E-04 

. 2851E-06 

. 1139E+02 

. 9868E+01 

.  1068E-04 

.  2  599E-06 

. 1042E+02 

. 8681E+01 

.  1066E-04 

.  2385E-06 

• 

. 9609E+01 

. 7650E+01 

.  1 065E-04 

.  2200E-06 

v 

.891 3E+01 

. 6739E+01 

. 1064E-04 

.  2038E-06 

V  .  * 

.831 IE+01 

. 5924E+01 

. 1063E-04 

.  1894E-06 

X 

. 7785E+01 

. 5 1 83E+0 1 

.  1062E-04 

. 1765E-06 

•  V- 

.  7321E+01 

. 4503E+01 

.  106  IE-04 

.  1 649E-06 

— -- 

.  6910E+01 

. 3870E+01 

.  1060E-04 

. 1545E-06 

.  6542E+01 

. 3276E+01 

. 1060E-04 

. 1450E-06 

.621 IE+01 

. 27 14E+01 

.  1059E-04 

. 1364E-06 

. 5912E+01 

.217  5E+0 1 

. 1058E-04 

.  1286E-06 

. 5641E+01 

.  1655E+01 

.  1058E-04 

.  1217E-06 

. 5393E+01 

.  1 149E+01 

. 1057E-04 

. 1157E-06 

w*  •  '  *  ' 

. 5166E+01 

. 6525E+00 

. 1056E-04 

.  1 105E-06 

w.  .  , 

. 4958E+01 

.  1612E+00 

. 1056E-04 

. 1063E-06 

. 4766E+01 

- . 3289E+00 

. 1055E-04 

. 1031E-06 

'v**'*-’* 

. 4587E+01 

-  .8162E+00 

. 1055E-04 

.  1011E-06 

»_*  »#*  •/ 

.  4422E+01 

- . 1297E+01 

. 1054E-04 

.  1001E-06 

.  4268E+01 

-.  1768E+01 

. 1054E-04 

.  1003E-06 

.  4 1 25E+01 

-  .  2226E+01 

. 1053E-04 

. 1013E-06 

-7 

. 3991E+01 

- . 2667E+01 

. 1053E-04 

. 1031E-06 

. 3865E+01 

-  .  3089E+01 

. 1052E-04 

. 1055E-06 

.  3747E+01 

- .  3491E+01 

. 1052E-04 

. 1083E-06 

.  3636E+01 

- .  3871E+01 

. 1051E-04 

. 1113E-06 

.  3531E+01 

- .  4227E+01 

. 1051E-04 

. 1 145E-06 

.  3432E+01 

- . 4560E+01 

. 1050E-04 

. 1 176E-06 

k-. 

.  3339E+01 

-  .  4868E+01 

. 1049E-04 

. 1206E-06 

/’v  V 

.  3250E+01 

- . 5154E+01 

. 1049E-04 

. 1236E-06 

.  3166E+01 

- .  541 6E+01 

. 1048E-04 

. 1263E-06 

.  3087E+01 

- . 5656E+01 

. 1048E-04 

. 1 289E-06 

y  *  »  "  l  ' 

.  301 IE+01 

-  .  5874E+01 

. 1047E-04 

. 1312E-06 

.  2939E+01 

- . 6073E+01 

. 1046E-04 

. 1334E-06 

. 1353E-06 
. 1371E-06 
. 1387E-06 
. 1400E-06 
. 1413E-06 
.  1424E-06 
. 1433E-06 
. 1441E-06 
. 1449E-06 
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Dispersion  Model  Input  Parameters 


Wind  Speed  =  5.0  m/sec 

Air  Temperature  =  293  °C 
Source  Rate  =  .38  kg/sec 
Stability  Class  =  F 
Liquid  Pool  Radius  =  9.8  m 


DOWNWIND  MAXIMUM  CLOUD  CLOUD  WIDTH  CLOUD  WIDTH  CLOUD  HEIGHT 

DISTANCE  CONCENTRATION  LOWER  FLAMMABLE  STOICHIOMETRIC  STOICHIOMETRIC 

LIMIT  MIXING  RATIO  MIXING  RATIO  MIXING  RATIO 

<  1.37.)  (  4.  Oy.)  <  4.07.) 

(METERS)  (V/V)  (METERS)  (METERS)  (METERS) 
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